home *** CD-ROM | disk | FTP | other *** search
/ System Booster / System Booster.iso / Archives / GNU / GNUPLOTsrc.lha / contour.c < prev    next >
C/C++ Source or Header  |  1996-01-22  |  45KB  |  1,287 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: contour.c,v 1.18 1995/12/02 21:16:34 drd Exp $";
  3. #endif
  4.  
  5.  
  6. /* GNUPLOT - contour.c */
  7. /*
  8.  * Copyright (C) 1986 - 1993   Thomas Williams, Colin Kelley
  9.  *
  10.  * Permission to use, copy, and distribute this software and its
  11.  * documentation for any purpose with or without fee is hereby granted, 
  12.  * provided that the above copyright notice appear in all copies and 
  13.  * that both that copyright notice and this permission notice appear 
  14.  * in supporting documentation.
  15.  *
  16.  * Permission to modify the software is granted, but not the right to
  17.  * distribute the modified code.  Modifications are to be distributed 
  18.  * as patches to released version.
  19.  *  
  20.  * This software is provided "as is" without express or implied warranty.
  21.  * 
  22.  *
  23.  * AUTHORS
  24.  * 
  25.  *   Original Software:
  26.  *       Gershon Elber
  27.  * 
  28.  * There is a mailing list for gnuplot users. Note, however, that the
  29.  * newsgroup 
  30.  *    comp.graphics.gnuplot 
  31.  * is identical to the mailing list (they
  32.  * both carry the same set of messages). We prefer that you read the
  33.  * messages through that newsgroup, to subscribing to the mailing list.
  34.  * (If you can read that newsgroup, and are already on the mailing list,
  35.  * please send a message info-gnuplot-request@dartmouth.edu, asking to be
  36.  * removed from the mailing list.)
  37.  *
  38.  * The address for mailing to list members is
  39.  *       info-gnuplot@dartmouth.edu
  40.  * and for mailing administrative requests is 
  41.  *       info-gnuplot-request@dartmouth.edu
  42.  * The mailing list for bug reports is 
  43.  *       bug-gnuplot@dartmouth.edu
  44.  * The list of those interested in beta-test versions is
  45.  *       info-gnuplot-beta@dartmouth.edu
  46.  */
  47.  
  48. #include <math.h> /* get prototype for sqrt */
  49. #include "plot.h"
  50. #include "setshow.h"
  51. #include <stdio.h>
  52.  
  53. #define DEFAULT_NUM_OF_ZLEVELS  10  /* Some dflt values (setable via flags). */
  54. #define DEFAULT_NUM_APPROX_PTS  5
  55. #define DEFAULT_BSPLINE_ORDER  3
  56. #define MAX_NUM_OF_ZLEVELS      100 /* Some max. values (setable via flags). */
  57. #define MAX_NUM_APPROX_PTS      100
  58. #define MAX_BSPLINE_ORDER      10
  59.  
  60. /* for some reason these symbols are also defined in plot.h under different */
  61. /* names */
  62. #define INTERP_NOTHING    CONTOUR_KIND_LINEAR    /* Kind of interpolations on contours. */
  63. #define INTERP_CUBIC    CONTOUR_KIND_CUBIC_SPL    /* Cubic spline interp. */
  64. #define APPROX_BSPLINE    CONTOUR_KIND_BSPLINE    /* Bspline interpolation. */
  65.  
  66. #define ACTIVE   1                    /* Status of edges at certain Z level. */
  67. #define INACTIVE 2
  68.  
  69. #define OPEN_CONTOUR     1                                 /* Contour kinds. */
  70. #define CLOSED_CONTOUR   2
  71.  
  72. #define EPSILON  1e-5              /* Used to decide if two float are equal. */
  73. #define INFINITY 1e10
  74.  
  75. #ifndef TRUE
  76. #define TRUE     -1
  77. #define FALSE    0
  78. #endif
  79.  
  80. #define DEFAULT_NUM_CONTOURS    10
  81. #define MAX_POINTS_PER_CNTR     100
  82. #define SHIFT_Z_EPSILON        0.000301060 /* Dec. change of poly bndry hit.*/
  83.  
  84. #define ABS(x)  ((x) > 0 ? (x) : (-(x)))
  85. #define SQR(x)  ((x) * (x))
  86.  
  87. typedef double table_entry[4];           /* Cubic spline interpolation 4 coef. */
  88.  
  89. struct vrtx_struct {
  90.     double X, Y, Z;                       /* The coordinates of this vertex. */
  91.     struct vrtx_struct *next;                             /* To chain lists. */
  92. };
  93.  
  94. struct edge_struct {
  95.     struct poly_struct *poly[2];   /* Each edge belongs to up to 2 polygons. */
  96.     struct vrtx_struct *vertex[2]; /* The two extreme points of this vertex. */
  97.     struct edge_struct *next;                             /* To chain lists. */
  98.     int status, /* Status flag to mark edges in scanning at certain Z level. */
  99.     boundary;                   /* True if this edge is on the boundary. */
  100. };
  101.  
  102. struct poly_struct {
  103.     struct edge_struct *edge[3];           /* As we do triangolation here... */
  104.     struct poly_struct *next;                             /* To chain lists. */
  105. };
  106.  
  107. struct cntr_struct {           /* Contours are saved using this struct list. */
  108.     double X, Y;                          /* The coordinates of this vertex. */
  109.     struct cntr_struct *next;                             /* To chain lists. */
  110. };
  111.  
  112. static int test_boundary;    /* If TRUE look for contours on boundary first. */
  113. static struct gnuplot_contours *contour_list = NULL;
  114. static double crnt_cntr[MAX_POINTS_PER_CNTR * 2];
  115. static int crnt_cntr_pt_index = 0;
  116. static double contour_level = 0.0;
  117. static table_entry *hermit_table = NULL;    /* Hold hermite table constants. */
  118. static int num_of_z_levels = DEFAULT_NUM_OF_ZLEVELS;  /* # Z contour levels. */
  119. static int num_approx_pts = DEFAULT_NUM_APPROX_PTS;/* # pts per approx/inter.*/
  120. static int bspline_order = DEFAULT_BSPLINE_ORDER;   /* Bspline order to use. */
  121. static int interp_kind = INTERP_NOTHING;  /* Linear, Cubic interp., Bspline. */
  122.  
  123.  
  124. static void add_cntr_point __P((double x, double y));
  125. static void end_crnt_cntr __P((void));
  126. static void gen_contours __P((struct edge_struct *p_edges, double z_level, double x_min, double x_max, double y_min, double y_max));
  127. static int update_all_edges __P((struct edge_struct *p_edges, double z_level));
  128. static struct cntr_struct *gen_one_contour __P((struct edge_struct *p_edges, double z_level, int *contour_kind, int *num_active));
  129. static struct cntr_struct *trace_contour __P((struct edge_struct *pe_start, double z_level, int *num_active, int contour_kind));
  130. static struct cntr_struct *update_cntr_pt __P((struct edge_struct *p_edge, double z_level, int *in_middle));
  131. static int fuzzy_equal __P((double x, double y));
  132. static void gen_triangle __P((int num_isolines, struct iso_curve *iso_lines,
  133. struct poly_struct **p_polys, struct edge_struct **p_edges, struct vrtx_struct
  134. **p_vrts, double *x_min, double *y_min, double *z_min, double *x_max, double
  135. *y_max, double *z_max));
  136. static struct vrtx_struct *gen_vertices __P((int grid_x_max, struct coordinate GPHUGE *points, double *x_min, double *y_min, double *z_min, double *x_max, double *y_max, double *z_max));
  137. static struct edge_struct *gen_edges __P((int grid_x_max, struct vrtx_struct *p_vrtx, struct edge_struct **pe_tail));
  138. static struct edge_struct *gen_edges_middle __P((int grid_x_max, struct vrtx_struct *p_vrtx1, struct vrtx_struct *p_vrtx2, struct edge_struct **pe_tail));
  139. static struct poly_struct *gen_polys __P((int grid_x_max, struct edge_struct *p_edge1, struct edge_struct *p_edge_middle, struct edge_struct *p_edge2, struct poly_struct **pp_tail));
  140. static void free_contour __P((struct cntr_struct *p_cntr));
  141. static void put_contour __P((struct cntr_struct *p_cntr, double z_level, double x_min, double x_max, double y_min, double y_max, int contr_kind));
  142. static void put_contour_nothing __P((struct cntr_struct *p_cntr));
  143. static void put_contour_cubic __P((struct cntr_struct *p_cntr, double z_level, double x_min, double x_max, double y_min, double y_max, int contr_kind));
  144. static void put_contour_bspline __P((struct cntr_struct *p_cntr, double z_level, double x_min, double x_max, double y_min, double y_max, int contr_kind));
  145. static void calc_tangent __P((int n, double x1, double x2, double x3, double y1, double y2, double y3, double *tx, double *ty));
  146. static int count_contour __P((struct cntr_struct *p_cntr));
  147. static void complete_spline_interp __P((struct cntr_struct *p_cntr, int n, double t_min, double t_max, double tx1, double ty1, double tx2, double ty2));
  148. static void calc_hermit_table __P((void));
  149. static void hermit_interp __P((double x1, double y1, double tx1, double ty1, double x2, double y2, double tx2, double ty2, double dt));
  150. static void prepare_spline_interp __P((double tangents_x[], double tangents_y[], struct cntr_struct *p_cntr, int n, double t_min, double t_max, double tx1, double ty1, double tx2, double ty2));
  151. static void gen_bspline_approx __P((struct cntr_struct *p_cntr, int num_of_points, int order, int contour_kind));
  152. static double fetch_knot __P((int contour_kind, int num_of_points, int order, int i));
  153. static void eval_bspline __P((double t, struct cntr_struct *p_cntr, int num_of_points, int order, int j, int contour_kind, double *x, double *y));
  154.  
  155. /*
  156.  * Entry routine to this whole set of contouring module.
  157.  */
  158. struct gnuplot_contours *contour(num_isolines, iso_lines,
  159.                  ZLevels, approx_pts, int_kind, order1,
  160.                  levels_kind, levels_list)
  161. int num_isolines;
  162. struct iso_curve *iso_lines;
  163. int ZLevels, approx_pts, int_kind, order1;
  164. int levels_kind;
  165. double *levels_list;
  166. {
  167.     int i;
  168.     struct poly_struct *p_polys, *p_poly;
  169.     struct edge_struct *p_edges, *p_edge;
  170.     struct vrtx_struct *p_vrts, *p_vrtx;
  171.     double x_min, y_min, z_min, x_max, y_max, z_max, z, dz;
  172.      struct gnuplot_contours *save_contour_list;
  173.  
  174.     num_of_z_levels = ZLevels;
  175.     num_approx_pts = approx_pts;
  176.     bspline_order = order1 - 1;
  177.     interp_kind = int_kind;
  178.  
  179.     contour_list = NULL;
  180.  
  181.     if (interp_kind == INTERP_CUBIC) calc_hermit_table();
  182.  
  183.     gen_triangle(num_isolines, iso_lines, &p_polys, &p_edges, &p_vrts,
  184.         &x_min, &y_min, &z_min, &x_max, &y_max, &z_max);
  185.     crnt_cntr_pt_index = 0;
  186.  
  187.     dz = (z_max - z_min) / (num_of_z_levels + 1);
  188.     z = z_min;
  189.     for (i = 0; i < num_of_z_levels; i++) {
  190.         switch(levels_kind) {
  191.             case LEVELS_AUTO:
  192.                 /* Step from z_min+dz upto z_max-dz in num_of_z_levels times. */
  193.                 z += dz;
  194.                 break;
  195.             case LEVELS_INCREMENTAL:
  196.                 z = levels_list[0] + i * levels_list[1];
  197.                 break;
  198.             case LEVELS_DISCRETE:
  199.                 z = is_log_z ? log(levels_list[i])/log_base_log_z : levels_list[i];
  200.                 break;
  201.         }
  202.         contour_level = z;
  203.          save_contour_list = contour_list;
  204.         gen_contours(p_edges, z + dz * SHIFT_Z_EPSILON, x_min, x_max,
  205.                                 y_min, y_max);
  206.          if(contour_list != save_contour_list) {
  207.              contour_list->isNewLevel = 1;
  208.              sprintf(contour_list->label, contour_format, is_log_z ? pow(base_log_z, z) : z);
  209.          }
  210.     }
  211.  
  212.     /* Free all contouring related temporary data. */
  213.     while (p_polys) {
  214.     p_poly = p_polys -> next;
  215.     free (p_polys);
  216.     p_polys = p_poly;
  217.     }
  218.     while (p_edges) {
  219.     p_edge = p_edges -> next;
  220.     free (p_edges);
  221.     p_edges = p_edge;
  222.     }
  223.     while (p_vrts) {
  224.     p_vrtx = p_vrts -> next;
  225.     free (p_vrts);
  226.     p_vrts = p_vrtx;
  227.     }
  228.  
  229.     if (interp_kind == INTERP_CUBIC) free(hermit_table);
  230.  
  231.     return contour_list;
  232. }
  233.  
  234. /*
  235.  * Adds another point to the currently build contour.
  236.  */
  237. static void add_cntr_point(x, y)
  238. double x, y;
  239. {
  240.     int index;
  241.  
  242.     if (crnt_cntr_pt_index >= MAX_POINTS_PER_CNTR-1) {
  243.     index = crnt_cntr_pt_index - 1;
  244.     end_crnt_cntr();
  245.     crnt_cntr[0] = crnt_cntr[index * 2];
  246.     crnt_cntr[1] = crnt_cntr[index * 2 + 1];
  247.     crnt_cntr_pt_index = 1; /* Keep the last point as first of this one. */
  248.     }
  249.     crnt_cntr[crnt_cntr_pt_index * 2] = x;
  250.     crnt_cntr[crnt_cntr_pt_index * 2 + 1] = y;
  251.     crnt_cntr_pt_index++;
  252. }
  253.  
  254. /*
  255.  * Done with current contour - create gnuplot data structure for it.
  256.  */
  257. static void end_crnt_cntr()
  258. {
  259.     int i;
  260.     struct gnuplot_contours *cntr = (struct gnuplot_contours *)
  261.                     alloc((unsigned long)sizeof(struct gnuplot_contours),
  262.                           "gnuplot_contour");
  263.  
  264.     cntr->coords = (struct coordinate GPHUGE *) alloc((unsigned long)sizeof(struct coordinate) *
  265.                               ( unsigned long)crnt_cntr_pt_index,
  266.                            "contour coords");
  267.     for (i=0; i<crnt_cntr_pt_index; i++) {
  268.     cntr->coords[i].x = crnt_cntr[i * 2];
  269.     cntr->coords[i].y = crnt_cntr[i * 2 + 1];
  270.     cntr->coords[i].z = contour_level;
  271.     }
  272.     cntr->num_pts = crnt_cntr_pt_index;
  273.  
  274.     cntr->next = contour_list;
  275.     contour_list = cntr;
  276.      contour_list->isNewLevel = 0;
  277.  
  278.     crnt_cntr_pt_index = 0;
  279. }
  280.  
  281. /*
  282.  * Generates all contours by tracing the intersecting triangles.
  283.  */
  284. static void gen_contours(p_edges, z_level, x_min, x_max, y_min, y_max)
  285. struct edge_struct *p_edges;
  286. double z_level, x_min, x_max, y_min, y_max;
  287. {
  288.     int num_active,                        /* Number of edges marked ACTIVE. */
  289.     contour_kind;                /* One of OPEN_CONTOUR, CLOSED_CONTOUR. */
  290.     struct cntr_struct *p_cntr;
  291.  
  292.     num_active = update_all_edges(p_edges, z_level);           /* Do pass 1. */
  293.  
  294.     test_boundary = TRUE;        /* Start to look for contour on boundaries. */
  295.  
  296.     while (num_active > 0) {                                   /* Do Pass 2. */
  297.         /* Generate One contour (and update MumActive as needed): */
  298.     p_cntr = gen_one_contour(p_edges, z_level, &contour_kind, &num_active);
  299.     put_contour(p_cntr, z_level, x_min, x_max, y_min, y_max,
  300.                   contour_kind); /* Emit it in requested format. */
  301.     }
  302. }
  303.  
  304. /*
  305.  * Does pass 1, or marks the edges which are active (crosses this z_level)
  306.  * as ACTIVE, and the others as INACTIVE:
  307.  * Returns number of active edges (marked ACTIVE).
  308.  */
  309. static int update_all_edges(p_edges, z_level)
  310. struct edge_struct *p_edges;
  311. double z_level;
  312. {
  313.     int count = 0;
  314.  
  315.     while (p_edges) {
  316.     if (((p_edges -> vertex[0] -> Z >= z_level) &&
  317.          (p_edges -> vertex[1] -> Z <= z_level)) ||
  318.         ((p_edges -> vertex[1] -> Z >= z_level) &&
  319.          (p_edges -> vertex[0] -> Z <= z_level))) {
  320.         p_edges -> status = ACTIVE;
  321.         count++;
  322.     }
  323.     else p_edges -> status = INACTIVE;
  324.     p_edges = p_edges -> next;
  325.     }
  326.  
  327.     return count;
  328. }
  329.  
  330. /*
  331.  * Does pass 2, or find one complete contour out of the triangolation data base:
  332.  * Returns a pointer to the contour (as linked list), contour_kind is set to
  333.  * one of OPEN_CONTOUR, CLOSED_CONTOUR, and num_active is updated.
  334.  */
  335. static struct cntr_struct *gen_one_contour(p_edges, z_level, contour_kind,
  336.                                 num_active)
  337. struct edge_struct *p_edges;
  338. double z_level;
  339. int *contour_kind, *num_active;
  340. {
  341.     struct edge_struct *pe_temp;
  342.  
  343.     if (test_boundary) {    /* Look for something to start with on boundary: */
  344.     pe_temp = p_edges;
  345.     while (pe_temp) {
  346.         if ((pe_temp -> status == ACTIVE) && (pe_temp -> boundary)) break;
  347.         pe_temp = pe_temp -> next;
  348.     }
  349.     if (!pe_temp) test_boundary = FALSE;/* No more contours on boundary. */
  350.     else {
  351.         *contour_kind = OPEN_CONTOUR;
  352.         return trace_contour(pe_temp, z_level, num_active, *contour_kind);
  353.     }
  354.     }
  355.  
  356.     if (!test_boundary) {        /* Look for something to start with inside: */
  357.     pe_temp = p_edges;
  358.     while (pe_temp) {
  359.         if ((pe_temp -> status == ACTIVE) && (!(pe_temp -> boundary)))
  360.         break;
  361.         pe_temp = pe_temp -> next;
  362.     }
  363.     if (!pe_temp) {
  364.         *num_active = 0;
  365.         return NULL;
  366.     }
  367.     else {
  368.         *contour_kind = CLOSED_CONTOUR;
  369.         return trace_contour(pe_temp, z_level, num_active, *contour_kind);
  370.     }
  371.     }
  372.     return NULL;             /* We should never be here, but lint... */
  373. }
  374.  
  375. /*
  376.  * Search the data base along a contour starts at the edge pe_start until
  377.  * a boundary edge is detected or until we close the loop back to pe_start.
  378.  * Returns a linked list of all the points on the contour
  379.  * Also decreases num_active by the number of points on contour.
  380.  */
  381. static struct cntr_struct *trace_contour(pe_start, z_level, num_active,
  382.                                 contour_kind)
  383. struct edge_struct *pe_start;
  384. double z_level;
  385. int *num_active, contour_kind;
  386. {
  387.     int i, in_middle;       /* If TRUE the z_level is in the middle of edge. */
  388.     struct cntr_struct *p_cntr, *pc_tail;
  389.     struct edge_struct *p_edge = pe_start, *p_next_edge;
  390.     struct poly_struct *p_poly, *PLastpoly = NULL;
  391.  
  392.     /* Generate the header of the contour - the point on pe_start. */
  393.     if (contour_kind == OPEN_CONTOUR) pe_start -> status = INACTIVE;
  394.     (*num_active)--;
  395.     p_cntr = pc_tail = update_cntr_pt(pe_start, z_level, &in_middle);
  396.     if (!in_middle) {
  397.     return NULL;
  398.     }
  399.  
  400.     do {
  401.     /* Find polygon to continue (Not where we came from - PLastpoly): */
  402.     if (p_edge -> poly[0] == PLastpoly) p_poly = p_edge -> poly[1];
  403.     else p_poly = p_edge -> poly[0];
  404.     p_next_edge = NULL;          /* In case of error, remains NULL. */
  405.     for (i=0; i<3; i++)              /* Test the 3 edges of the polygon: */
  406.         if (p_poly -> edge[i] != p_edge)
  407.             if (p_poly -> edge[i] -> status == ACTIVE)
  408.             p_next_edge = p_poly -> edge[i];
  409.     if (!p_next_edge) {
  410.         pc_tail -> next = NULL;
  411.         free_contour(p_cntr);
  412.         return NULL;
  413.     }
  414.     p_edge = p_next_edge;
  415.     PLastpoly = p_poly;
  416.     p_edge -> status = INACTIVE;
  417.     (*num_active)--;
  418.     pc_tail -> next = update_cntr_pt(p_edge, z_level, &in_middle);
  419.     if (!in_middle) {
  420.         pc_tail -> next = NULL;
  421.         free_contour(p_cntr);
  422.         return NULL;
  423.     }
  424.         pc_tail = pc_tail -> next;
  425.     }
  426.     while ((pe_start != p_edge) && (!p_edge -> boundary));
  427.     pc_tail -> next = NULL;
  428.  
  429.     return p_cntr;
  430. }
  431.  
  432. /*
  433.  * Allocates one contour location and update it to to correct position
  434.  * according to z_level and edge p_edge. if z_level is found to be at
  435.  * one of the extreme points nothing is allocated (NULL is returned)
  436.  * and in_middle is set to FALSE.
  437.  */
  438. static struct cntr_struct *update_cntr_pt(p_edge, z_level, in_middle)
  439. struct edge_struct *p_edge;
  440. double z_level;
  441. int *in_middle;
  442. {
  443.     double t;
  444.     struct cntr_struct *p_cntr;
  445.  
  446.     t = (z_level - p_edge -> vertex[0] -> Z) /
  447.     (p_edge -> vertex[1] -> Z - p_edge -> vertex[0] -> Z);
  448.  
  449.     if (fuzzy_equal(t, 1.0) || fuzzy_equal(t, 0.0)) {
  450.         *in_middle = FALSE;
  451.         return NULL;
  452.     }
  453.     else {
  454.     *in_middle = TRUE;
  455.     p_cntr = (struct cntr_struct *) alloc((unsigned long)sizeof(struct cntr_struct),
  456.                             "contour cntr_struct");
  457.     p_cntr -> X = p_edge -> vertex[1] -> X * t +
  458.              p_edge -> vertex[0] -> X * (1-t);
  459.     p_cntr -> Y = p_edge -> vertex[1] -> Y * t +
  460.              p_edge -> vertex[0] -> Y * (1-t);
  461.     return p_cntr;
  462.     }
  463. }
  464.  
  465. /*
  466.  * Simple routine to decide if two real values are equal by simply
  467.  * calculating the relative/absolute error between them (< EPSILON).
  468.  */
  469. static int fuzzy_equal(x, y)
  470. double x, y;
  471. {
  472.     if (ABS(x) > EPSILON)            /* Calculate relative error: */
  473.         return (ABS((x - y) / x) < EPSILON);
  474.     else                    /* Calculate absolute error: */
  475.     return (ABS(x - y) < EPSILON);
  476. }
  477.  
  478. /*
  479.  * Generate the triangles.
  480.  * Returns the lists (vrtxs edges & polys) via pointers to their heads.
  481.  */
  482. static void gen_triangle(num_isolines, iso_lines, p_polys, p_edges,
  483.     p_vrts, x_min, y_min, z_min, x_max, y_max, z_max)
  484. int num_isolines;
  485. struct iso_curve *iso_lines;
  486. struct poly_struct **p_polys;
  487. struct edge_struct **p_edges;
  488. struct vrtx_struct **p_vrts;
  489. double *x_min, *y_min, *z_min, *x_max, *y_max, *z_max;
  490. {
  491.     int i, grid_x_max = iso_lines->p_count;
  492.     struct vrtx_struct *p_vrtx1, *p_vrtx2, *pv_temp;
  493.     struct edge_struct *p_edge1, *p_edge2, *pe_tail1, *pe_tail2, *pe_temp,
  494.               *p_edge_middle, *pe_m_tail;
  495.     struct poly_struct *p_poly, *pp_tail;
  496.  
  497.     *p_polys = NULL;
  498.     *p_edges = NULL;
  499.     *p_vrts = NULL;
  500.     *z_min = INFINITY;
  501.     *y_min = INFINITY;
  502.     *x_min = INFINITY;
  503.     *z_max = -INFINITY;
  504.     *y_max = -INFINITY;
  505.     *x_max = -INFINITY;
  506.  
  507.     /* Read 1st row. */
  508.     p_vrtx1 = gen_vertices(grid_x_max, iso_lines->points,
  509.                x_min, y_min, z_min, x_max, y_max, z_max);
  510.     *p_vrts = p_vrtx1;
  511.     /* Gen. its edges.*/
  512.     pe_temp = p_edge1 = gen_edges(grid_x_max, p_vrtx1, &pe_tail1);
  513.     for (i = 1; i < grid_x_max; i++) {/* Mark one side of edges as boundary. */
  514.     pe_temp -> poly[1] = NULL;
  515.     pe_temp = pe_temp -> next;
  516.     }
  517.     for (i = 1; i < num_isolines; i++) { /* Read next column and gen. polys. */
  518.     iso_lines = iso_lines->next;
  519.         /* Get row into list. */
  520.         p_vrtx2 = gen_vertices(grid_x_max, iso_lines->points,
  521.                    x_min, y_min, z_min, x_max, y_max, z_max);
  522.         /* Generate its edges. */
  523.         p_edge2 = gen_edges(grid_x_max, p_vrtx2, &pe_tail2);
  524.     /* Generate edges from one vertex list to the other one: */
  525.     p_edge_middle = gen_edges_middle(grid_x_max, p_vrtx1, p_vrtx2,
  526.                                  &pe_m_tail);
  527.  
  528.     /* Now we can generate the polygons themselves (triangles). */
  529.     p_poly = gen_polys(grid_x_max, p_edge1, p_edge_middle, p_edge2,
  530.                                  &pp_tail);
  531.         pe_tail1 -> next = (*p_edges);      /* Chain new edges to main list. */
  532.         pe_m_tail -> next = p_edge1;
  533.     *p_edges = p_edge_middle;
  534.     pe_tail1 = pe_tail2;
  535.     p_edge1 = p_edge2;
  536.  
  537.     pv_temp = p_vrtx2;
  538.     while (pv_temp -> next) pv_temp = pv_temp -> next;
  539.     pv_temp -> next = *p_vrts;
  540.     *p_vrts = p_vrtx1 = p_vrtx2;
  541.  
  542.         pp_tail -> next = (*p_polys);       /* Chain new polys to main list. */
  543.     *p_polys = p_poly;
  544.     }
  545.  
  546.     pe_temp = p_edge1;
  547.     for (i = 1; i < grid_x_max; i++) {/* Mark one side of edges as boundary. */
  548.     pe_temp -> poly[0] = NULL;
  549.     pe_temp = pe_temp -> next;
  550.     }
  551.  
  552.     pe_tail1 -> next = (*p_edges);    /* Chain last edges list to main list. */
  553.     *p_edges = p_edge1;
  554.  
  555.     /* Update the boundary flag, saved in each edge, and update indexes: */
  556.     pe_temp = (*p_edges);
  557.     i = 1;
  558.  
  559.     while (pe_temp) {
  560.     pe_temp -> boundary = (!(pe_temp -> poly[0])) ||
  561.                   (!(pe_temp -> poly[1]));
  562.     pe_temp = pe_temp -> next;
  563.     }
  564. }
  565.  
  566. /*
  567.  * Handles grid_x_max 3D points (One row) and generate linked list for them.
  568.  */
  569. static struct vrtx_struct *gen_vertices(grid_x_max, points,
  570.                       x_min, y_min, z_min, x_max, y_max, z_max)
  571. int grid_x_max;
  572. struct coordinate GPHUGE *points;
  573. double *x_min, *y_min, *z_min, *x_max, *y_max, *z_max;
  574. {
  575.     int i;
  576.      /* avoid gcc -Wall warning about unused variables */
  577.     struct vrtx_struct *p_vrtx=NULL, *pv_tail=NULL, *pv_temp;
  578.  
  579.     for (i=0; i<grid_x_max; i++) {/* Get a point and generate the structure. */
  580.         pv_temp = (struct vrtx_struct *) alloc((unsigned long)sizeof(struct vrtx_struct),
  581.                         "contour vertex");
  582.     pv_temp -> X = points[i].x;
  583.     pv_temp -> Y = points[i].y;
  584.     pv_temp -> Z = points[i].z;
  585.  
  586.     if (pv_temp -> X > *x_max) *x_max = pv_temp -> X; /* Update min/max. */
  587.     if (pv_temp -> Y > *y_max) *y_max = pv_temp -> Y;
  588.     if (pv_temp -> Z > *z_max) *z_max = pv_temp -> Z;
  589.     if (pv_temp -> X < *x_min) *x_min = pv_temp -> X;
  590.     if (pv_temp -> Y < *y_min) *y_min = pv_temp -> Y;
  591.     if (pv_temp -> Z < *z_min) *z_min = pv_temp -> Z;
  592.  
  593.     if (i == 0)                              /* First vertex in row: */
  594.         p_vrtx = pv_tail = pv_temp;
  595.     else {
  596.         pv_tail -> next = pv_temp;   /* Stick new record as last one. */
  597.         pv_tail = pv_tail -> next;    /* And continue to last record. */
  598.     }
  599.     }
  600.     pv_tail -> next = NULL;
  601.  
  602.     return p_vrtx;
  603. }
  604.  
  605. /*
  606.  * Combines N vertices in pair to form N-1 edges.
  607.  * Returns pointer to the edge list (pe_tail will point on last edge in list).
  608.  */
  609. static struct edge_struct *gen_edges(grid_x_max, p_vrtx, pe_tail)
  610. int grid_x_max;
  611. struct vrtx_struct *p_vrtx;
  612. struct edge_struct **pe_tail;
  613. {
  614.     int i;
  615.     struct edge_struct *p_edge=NULL, *pe_temp;
  616.  
  617.     for (i=0; i<grid_x_max-1; i++) {         /* Generate grid_x_max-1 edges: */
  618.     pe_temp = (struct edge_struct *) alloc((unsigned long)sizeof(struct edge_struct),
  619.                         "contour edge");
  620.     pe_temp -> vertex[0] = p_vrtx;              /* First vertex of edge. */
  621.     p_vrtx = p_vrtx -> next;                     /* Skip to next vertex. */
  622.     pe_temp -> vertex[1] = p_vrtx;             /* Second vertex of edge. */
  623.         if (i == 0)                                    /* First edge in row: */
  624.         p_edge = (*pe_tail) = pe_temp;
  625.     else {
  626.         (*pe_tail) -> next = pe_temp;   /* Stick new record as last one. */
  627.         *pe_tail = (*pe_tail) -> next;   /* And continue to last record. */
  628.      }
  629.     }
  630.     (*pe_tail) -> next = NULL;
  631.  
  632.     return p_edge;
  633. }
  634.  
  635. /*
  636.  * Combines 2 lists of N vertices each into edge list:
  637.  * The dots (.) are the vertices list, and the              .  .  .  .
  638.  *  edges generated are alternations of vertical edges      |\ |\ |\ |
  639.  *  (|) and diagonal ones (\).                              | \| \| \|
  640.  *  A pointer to edge list (alternate | , \) is returned    .  .  .  .
  641.  * Note this list will have (2*grid_x_max-1) edges (pe_tail points on last
  642.  * record).
  643.  */
  644. static struct edge_struct *gen_edges_middle(grid_x_max, p_vrtx1, p_vrtx2,
  645.                                 pe_tail)
  646. int grid_x_max;
  647. struct vrtx_struct *p_vrtx1, *p_vrtx2;
  648. struct edge_struct **pe_tail;
  649. {
  650.     int i;
  651.     struct edge_struct *p_edge, *pe_temp;
  652.  
  653.     /* Gen first (|). */
  654.     pe_temp = (struct edge_struct *) alloc((unsigned long)sizeof(struct edge_struct),
  655.                             "contour edge");
  656.     pe_temp -> vertex[0] = p_vrtx2;                 /* First vertex of edge. */
  657.     pe_temp -> vertex[1] = p_vrtx1;                /* Second vertex of edge. */
  658.     p_edge = (*pe_tail) = pe_temp;
  659.  
  660.     /* Advance in vrtx list grid_x_max-1 times, and gen. 2 edges /| for each.*/
  661.     for (i=0; i<grid_x_max-1; i++) {
  662.     /* The / edge. */
  663.     pe_temp = (struct edge_struct *) alloc((unsigned long)sizeof(struct edge_struct),
  664.                             "contour edge");
  665.     pe_temp -> vertex[0] = p_vrtx1;             /* First vertex of edge. */
  666.     pe_temp -> vertex[1] = p_vrtx2 -> next;    /* Second vertex of edge. */
  667.         (*pe_tail) -> next = pe_temp;       /* Stick new record as last one. */
  668.     *pe_tail = (*pe_tail) -> next;       /* And continue to last record. */
  669.  
  670.     /* The | edge. */
  671.     pe_temp = (struct edge_struct *) alloc((unsigned long)sizeof(struct edge_struct),
  672.                             "contour edge");
  673.     pe_temp -> vertex[0] = p_vrtx2 -> next;     /* First vertex of edge. */
  674.     pe_temp -> vertex[1] = p_vrtx1 -> next;    /* Second vertex of edge. */
  675.         (*pe_tail) -> next = pe_temp;       /* Stick new record as last one. */
  676.     *pe_tail = (*pe_tail) -> next;       /* And continue to last record. */
  677.  
  678.         p_vrtx1 = p_vrtx1 -> next;   /* Skip to next vertices in both lists. */
  679.         p_vrtx2 = p_vrtx2 -> next;
  680.     }
  681.     (*pe_tail) -> next = NULL;
  682.  
  683.     return p_edge;
  684. }
  685.  
  686. /*
  687.  * Combines 3 lists of edges into triangles:
  688.  * 1. p_edge1: Top horizontal edge list:        -----------------------
  689.  * 2. p_edge_middge: middle edge list:         |\  |\  |\  |\  |\  |\  |
  690.  *                                             |  \|  \|  \|  \|  \|  \|
  691.  * 3. p_edge2: Bottom horizontal edge list:     -----------------------
  692.  * Note that p_edge1/2 lists has grid_x_max-1 edges, while p_edge_middle has
  693.  * (2*grid_x_max-1) edges.
  694.  * The routine simple scans the two list    Upper 1         Lower
  695.  * and generate two triangle upper one        ----         | \
  696.  * and lower one from the lists:             0\   |2      0|   \1
  697.  * (Nums. are edges order in polys)             \ |         ----
  698.  * The routine returns a pointer to a                         2
  699.  * polygon list (pp_tail points on last polygon).          1
  700.  *                                                   -----------
  701.  * In addition, the edge lists are updated -        | \   0     |
  702.  * each edge has two pointers on the two            |   \       |
  703.  * (one active if boundary) polygons which         0|1   0\1   0|1
  704.  * uses it. These two pointer to polygons           |       \   |
  705.  * are named: poly[0], poly[1]. The diagram         |    1    \ |
  706.  * on the right show how they are used for the       -----------
  707.  * upper and lower polygons.                             0
  708.  */
  709. static struct poly_struct *gen_polys(grid_x_max, p_edge1, p_edge_middle,
  710.                             p_edge2, pp_tail)
  711. int grid_x_max;
  712. struct edge_struct *p_edge1, *p_edge_middle, *p_edge2;
  713. struct poly_struct **pp_tail;
  714. {
  715.     int i;
  716.     struct poly_struct *p_poly=NULL, *pp_temp;
  717.  
  718.     p_edge_middle -> poly[0] = NULL;                /* Its boundary! */
  719.  
  720.     /* Advance in vrtx list grid_x_max-1 times, and gen. 2 polys for each. */
  721.     for (i=0; i<grid_x_max-1; i++) {
  722.     /* The Upper. */
  723.     pp_temp = (struct poly_struct *) alloc((unsigned long)sizeof(struct poly_struct),
  724.                             "contour poly");
  725.     /* Now update polys about its edges, and edges about the polygon. */
  726.     pp_temp -> edge[0] = p_edge_middle -> next;
  727.     p_edge_middle -> next -> poly[1] = pp_temp;
  728.     pp_temp -> edge[1] = p_edge1;
  729.     p_edge1 -> poly[0] = pp_temp;
  730.     pp_temp -> edge[2] = p_edge_middle -> next -> next;
  731.     p_edge_middle -> next -> next -> poly[0] = pp_temp;
  732.     if (i == 0)                   /* Its first one in list: */
  733.         p_poly = (*pp_tail) = pp_temp;
  734.     else {
  735.         (*pp_tail) -> next = pp_temp;
  736.         *pp_tail = (*pp_tail) -> next;
  737.     }
  738.  
  739.     /* The Lower. */
  740.     pp_temp = (struct poly_struct *) alloc((unsigned long)sizeof(struct poly_struct),
  741.                             "contour poly");
  742.     /* Now update polys about its edges, and edges about the polygon. */
  743.     pp_temp -> edge[0] = p_edge_middle;
  744.     p_edge_middle -> poly[1] = pp_temp;
  745.     pp_temp -> edge[1] = p_edge_middle -> next;
  746.     p_edge_middle -> next -> poly[0] = pp_temp;
  747.     pp_temp -> edge[2] = p_edge2;
  748.     p_edge2 -> poly[1] = pp_temp;
  749.     (*pp_tail) -> next = pp_temp;
  750.     *pp_tail = (*pp_tail) -> next;
  751.  
  752.         p_edge1 = p_edge1 -> next;
  753.         p_edge2 = p_edge2 -> next;
  754.         p_edge_middle = p_edge_middle -> next -> next;
  755.     }
  756.     p_edge_middle -> poly[1] = NULL;                /* Its boundary! */
  757.     (*pp_tail) -> next = NULL;
  758.  
  759.     return p_poly;
  760. }
  761.  
  762. /*
  763.  * Calls the (hopefully) desired interpolation/approximation routine.
  764.  */
  765. static void put_contour(p_cntr, z_level, x_min, x_max, y_min, y_max, contr_kind)
  766. struct cntr_struct *p_cntr;
  767. double z_level, x_min, x_max, y_min, y_max;
  768. int contr_kind;
  769. {
  770.     if (!p_cntr) return;            /* Nothing to do if it is empty contour. */
  771.  
  772.     switch (interp_kind) {
  773.     case INTERP_NOTHING:              /* No interpolation/approximation. */
  774.         put_contour_nothing(p_cntr);
  775.         break;
  776.     case INTERP_CUBIC:                    /* Cubic spline interpolation. */
  777.         put_contour_cubic(p_cntr, z_level, x_min, x_max, y_min, y_max,
  778.                                 contr_kind);
  779.         break;
  780.     case APPROX_BSPLINE:                       /* Bspline approximation. */
  781.         put_contour_bspline(p_cntr, z_level, x_min, x_max, y_min, y_max,
  782.                                     contr_kind);
  783.         break;
  784.     }
  785.  
  786.     free_contour(p_cntr);
  787. }
  788.  
  789. /*
  790.  * Simply puts contour coordinates in order with no interpolation or
  791.  * approximation.
  792.  */
  793. static void put_contour_nothing(p_cntr)
  794. struct cntr_struct *p_cntr;
  795. {
  796.     while (p_cntr) {
  797.     add_cntr_point(p_cntr -> X, p_cntr -> Y);
  798.     p_cntr = p_cntr -> next;
  799.     }
  800.     end_crnt_cntr();
  801. }
  802.  
  803. /*
  804.  * Find Complete Cubic Spline Interpolation.
  805.  */
  806. static void put_contour_cubic(p_cntr, z_level, x_min, x_max, y_min, y_max,
  807.                                  contr_kind)
  808. struct cntr_struct *p_cntr;
  809. double z_level, x_min, x_max, y_min, y_max;
  810. int contr_kind;
  811. {
  812.     int num_pts, i;
  813.     double tx1, ty1, tx2, ty2;                    /* Tangents at end points. */
  814.     struct cntr_struct *pc_temp;
  815.  
  816.     num_pts = count_contour(p_cntr);         /* Number of points in contour. */
  817.  
  818.     if (num_pts > 2) {  /* Take into account 3 points in tangent estimation. */
  819.     calc_tangent(3, p_cntr -> X, p_cntr -> next -> X,
  820.             p_cntr -> next -> next -> X,
  821.             p_cntr -> Y, p_cntr -> next -> Y,
  822.             p_cntr -> next -> next -> Y, &tx1, &ty1);
  823.     pc_temp = p_cntr;
  824.     for (i=3; i<num_pts; i++) pc_temp = pc_temp -> next;/* Go to the end.*/
  825.     calc_tangent(3, pc_temp -> next -> next -> X,
  826.              pc_temp -> next -> X, pc_temp -> X,
  827.             pc_temp -> next -> next -> Y,
  828.             pc_temp -> next -> Y, pc_temp -> Y, &tx2, &ty2);
  829.         tx2 = (-tx2);   /* Inverse the vector as we need opposite direction. */
  830.         ty2 = (-ty2);
  831.     }
  832.     /* If following (num_pts > 1) is TRUE then exactly 2 points in contour.  */
  833.     else if (num_pts > 1) {/* Take into account 2 points in tangent estimat. */
  834.     calc_tangent(2, p_cntr -> X, p_cntr -> next -> X, 0.0,
  835.             p_cntr -> Y, p_cntr -> next -> Y, 0.0, &tx1, &ty1);
  836.     calc_tangent(2, p_cntr -> next -> X, p_cntr -> X, 0.0,
  837.             p_cntr -> next -> Y, p_cntr -> Y, 0.0, &tx2, &ty2);
  838.         tx2 = (-tx2);   /* Inverse the vector as we need opposite direction. */
  839.         ty2 = (-ty2);
  840.     }
  841.     else return;            /* Only one point ( ?? ) - ignore it. */
  842.  
  843.     switch (contr_kind) {
  844.     case OPEN_CONTOUR:
  845.         break;
  846.     case CLOSED_CONTOUR:
  847.         tx1 = tx2 = (tx1 + tx2) / 2.0;         /* Make tangents equal. */
  848.         ty1 = ty2 = (ty1 + ty2) / 2.0;
  849.         break;
  850.     }
  851.     complete_spline_interp(p_cntr, num_pts, 0.0, 1.0, tx1, ty1, tx2, ty2);
  852.     end_crnt_cntr();
  853. }
  854.  
  855. /*
  856.  * Find Bspline approximation for this data set.
  857.  * Uses global variable num_approx_pts to determine number of samples per
  858.  * interval, where the knot vector intervals are assumed to be uniform, and
  859.  * Global variable bspline_order for the order of Bspline to use.
  860.  */
  861. static void put_contour_bspline(p_cntr, z_level, x_min, x_max, y_min, y_max,
  862.                                 contr_kind)
  863. struct cntr_struct *p_cntr;
  864. double z_level, x_min, x_max,  y_min, y_max;
  865. int contr_kind;
  866. {
  867.     int num_pts, order = bspline_order;
  868.  
  869.     num_pts = count_contour(p_cntr);         /* Number of points in contour. */
  870.     if (num_pts < 2) return;     /* Can't do nothing if empty or one points! */
  871.     /* Order must be less than number of points in curve - fix it if needed. */
  872.     if (order > num_pts - 1) order = num_pts - 1;
  873.  
  874.     gen_bspline_approx(p_cntr, num_pts, order, contr_kind);
  875.     end_crnt_cntr();
  876. }
  877.  
  878. /*
  879.  * Estimate the tangents according to the n last points where n might be
  880.  * 2 or 3 (if 2 onlt x1, x2).
  881.  */
  882. static void calc_tangent(n, x1, x2, x3, y1, y2, y3, tx, ty)
  883. int n;
  884. double x1, x2, x3, y1, y2, y3, *tx, *ty;
  885. {
  886.     double v1[2], v2[2], v1_magnitude, v2_magnitude;
  887.  
  888.     switch (n) {
  889.     case 2:
  890.         *tx = (x2 - x1) * 0.3;
  891.         *ty = (y2 - y1) * 0.3;
  892.         break;
  893.     case 3:
  894.         v1[0] = x2 - x1;   v1[1] = y2 - y1;
  895.         v2[0] = x3 - x2;   v2[1] = y3 - y2;
  896.         v1_magnitude = sqrt(SQR(v1[0]) + SQR(v1[1]));
  897.         v2_magnitude = sqrt(SQR(v2[0]) + SQR(v2[1]));
  898.         *tx = (v1[0] / v1_magnitude) - (v2[0] / v2_magnitude) * 0.1;
  899.         *tx *= v1_magnitude * 0.1;  /* Make tangent less than magnitude. */
  900.         *ty = (v1[1] / v1_magnitude) - (v2[1] / v2_magnitude) * 0.1;
  901.         *ty *= v1_magnitude * 0.1;  /* Make tangent less than magnitude. */
  902.         break;
  903.     default:                       /* Should not happen! */
  904.         (*ty) = 0.1;
  905.         *tx = 0.1;
  906.         break;
  907.     }
  908. }
  909.  
  910. /*
  911.  * Free all elements in the contour list.
  912.  */
  913. static void free_contour(p_cntr)
  914. struct cntr_struct *p_cntr;
  915. {
  916.     struct cntr_struct *pc_temp;
  917.  
  918.     while (p_cntr) {
  919.     pc_temp = p_cntr;
  920.     p_cntr = p_cntr -> next;
  921.     free((char *) pc_temp);
  922.     }
  923. }
  924.  
  925. /*
  926.  * Counts number of points in contour.
  927.  */
  928. static int count_contour(p_cntr)
  929. struct cntr_struct *p_cntr;
  930. {
  931.     int count = 0;
  932.  
  933.     while (p_cntr) {
  934.     count++;
  935.     p_cntr = p_cntr -> next;
  936.     }
  937.     return count;
  938. }
  939.  
  940. /*
  941.  * Interpolate given point list (defined via p_cntr) using Complete
  942.  * Spline interpolation.
  943.  */
  944. static void complete_spline_interp(p_cntr, n, t_min, t_max, tx1, ty1, tx2, ty2)
  945. struct cntr_struct *p_cntr;
  946. int n;
  947. double t_min, t_max, tx1, ty1, tx2, ty2;
  948. {
  949.     double dt, *tangents_x, *tangents_y;
  950.     int i;
  951.  
  952.     tangents_x = (double *) alloc((unsigned long) (sizeof(double) * n),
  953.                         "contour c_s_intr");
  954.     tangents_y = (double *) alloc((unsigned long) (sizeof(double) * n),
  955.                         "contour c_s_intr");
  956.  
  957.     if (n > 1) prepare_spline_interp(tangents_x, tangents_y, p_cntr, n,
  958.                     t_min, t_max, tx1, ty1, tx2, ty2);
  959.     else {
  960.     free((char *) tangents_x);
  961.     free((char *) tangents_y);
  962.     return;
  963.     }
  964.  
  965.     dt = (t_max-t_min)/(n-1);
  966.  
  967.     add_cntr_point(p_cntr -> X, p_cntr -> Y);           /* First point. */
  968.  
  969.     for (i=0; i<n-1; i++) {
  970.         hermit_interp(p_cntr -> X, p_cntr -> Y,
  971.                      tangents_x[i], tangents_y[i],
  972.              p_cntr -> next -> X, p_cntr -> next -> Y,
  973.                      tangents_x[i+1], tangents_y[i+1], dt);
  974.  
  975.         p_cntr = p_cntr -> next;
  976.     }
  977.  
  978.     free((char *) tangents_x);
  979.     free((char *) tangents_y);
  980. }
  981.  
  982. /*
  983.  * Routine to calculate intermidiate value of the Hermit Blending function:
  984.  * This routine should be called only ONCE at the beginning of the program.
  985.  */
  986. static void calc_hermit_table()
  987. {
  988.     int i;
  989.     double t, dt;
  990.  
  991.     hermit_table = (table_entry *) alloc ((unsigned long) (sizeof(table_entry) *
  992.                         (num_approx_pts + 1)),
  993.                         "contour hermit table");
  994.     t = 0;
  995.     dt = 1.0/num_approx_pts;
  996.     for (i=0; i<=num_approx_pts; i++) {
  997.         hermit_table[i][0] = (t-1)*(t-1)*(2*t+1);             /* h00. */
  998.         hermit_table[i][1] = t*t*(-2*t+3);                 /* h10. */
  999.         hermit_table[i][2] = t*(t-1)*(t-1);                 /* h01. */
  1000.         hermit_table[i][3] = t*t*(t-1);                     /* h11. */
  1001.         t = t + dt;
  1002.     }
  1003. }
  1004.  
  1005. /*
  1006.  * Routine to generate an hermit interpolation between two points given as
  1007.  * two InterpStruct structures. Assume hermit_table is already calculated.
  1008.  * Currently the points generated are printed to stdout as two reals (X, Y).
  1009.  */
  1010. static void hermit_interp(x1, y1, tx1, ty1, x2, y2, tx2, ty2, dt)
  1011. double x1, y1, tx1, ty1, x2, y2, tx2, ty2, dt;
  1012. {
  1013.     int i;
  1014.     double x, y, vec_size, tang_size;
  1015.  
  1016.     tx1 *= dt;  ty1 *= dt; /* Normalize the tangents according to param. t.  */
  1017.     tx2 *= dt;  ty2 *= dt;
  1018.  
  1019.     /* Normalize the tangents so that their magnitude will be 1/3 of the     */
  1020.     /* segment length. This tumb rule guaranteed no cusps or loops!          */
  1021.     /* Note that this normalization keeps continuity to be G1 (but not C1).  */
  1022.     vec_size = sqrt(SQR(x1 - x2) + SQR(y2 - y1));
  1023.     tang_size = sqrt(SQR(tx1) + SQR(ty1));                  /* Normalize T1. */
  1024.     if (tang_size * 3 > vec_size) {
  1025.     tx1 *= vec_size / (tang_size * 3);
  1026.     ty1 *= vec_size / (tang_size * 3);
  1027.     }
  1028.     tang_size = sqrt(SQR(tx2) + SQR(ty2));                  /* Normalize T2. */
  1029.     if (tang_size * 3 > vec_size) {
  1030.     tx2 *= vec_size / (tang_size * 3);
  1031.     ty2 *= vec_size / (tang_size * 3);
  1032.     }
  1033.  
  1034.     for (i=1; i<=num_approx_pts; i++) {      /* Note we start from 1 - first */
  1035.         x = hermit_table[i][0] * x1 +       /* point is not printed as it is */
  1036.             hermit_table[i][1] * x2 +   /* redundent (last on last section). */
  1037.             hermit_table[i][2] * tx1 +
  1038.             hermit_table[i][3] * tx2;
  1039.         y = hermit_table[i][0] * y1 +
  1040.             hermit_table[i][1] * y2 +
  1041.             hermit_table[i][2] * ty1 +
  1042.             hermit_table[i][3] * ty2;
  1043.     add_cntr_point(x, y);
  1044.     }
  1045. }
  1046.  
  1047. /*
  1048.  * Routine to Set up the 3*N mat for solve_tri_diag routine used in the
  1049.  * Complete Spline Interpolation. Returns TRUE of calc O.K.
  1050.  * Gets the points list in p_cntr (Of length n) and with tangent vectors tx1,
  1051.  * ty1 at starting point and tx2, ty2 and end point.
  1052.  */
  1053. static void prepare_spline_interp(tangents_x, tangents_y, p_cntr, n, t_min, t_max,
  1054.                tx1, ty1, tx2, ty2)
  1055. double tangents_x[], tangents_y[];
  1056. struct cntr_struct *p_cntr;
  1057. int n;
  1058. double t_min, t_max, tx1, ty1, tx2, ty2;
  1059. {
  1060.     int i;
  1061.     double *r, t, dt;
  1062.     tri_diag *m;           /* The tri-diagonal matrix is saved here. */
  1063.     struct cntr_struct *p;
  1064.  
  1065.     m = (tri_diag *) alloc((unsigned long) (sizeof(tri_diag) * n),
  1066.                         "contour tri_diag");
  1067.     r = (double *) alloc((unsigned long) (sizeof(double) * n),
  1068.                         "contour tri_diag2");
  1069.     n--;
  1070.  
  1071.     p = p_cntr;
  1072.     m[0][0] = 0.0;    m[0][1] = 1.0;    m[0][2] = 0.0;
  1073.     m[n][0] = 0.0;    m[n][1] = 1.0;    m[n][2] = 0.0;
  1074.     r[0] = tx1;                           /* Set start tangent. */
  1075.     r[n] = tx2;                         /* Set end tangent. */
  1076.     t = t_min;
  1077.     dt = (t_max-t_min)/n;
  1078.     for (i=1; i<n; i++) {
  1079.        t = t + dt;
  1080.        m[i][0] = dt;
  1081.        m[i][2] = dt;
  1082.        m[i][1] = 2 * (m[i][0] + m[i][2]);
  1083.        r[i] = m[i][0] * ((p -> next -> X) - (p -> X)) / m[i][2]
  1084.             + m[i][2] * ((p -> next -> next -> X) - 
  1085.                          (p -> next -> X)) / m[i][0];
  1086.        r[i] *= 3.0;
  1087.        p = p -> next;
  1088.     }
  1089.  
  1090.     if (!solve_tri_diag(m, r, tangents_x, n+1)) { /* Find the X(t) tangents. */
  1091.         free((char *) m);
  1092.         free((char *) r);
  1093.     int_error("Cannt interpolate X using complete splines", NO_CARET);
  1094.     }
  1095.  
  1096.     p = p_cntr;
  1097.     m[0][0] = 0.0;    m[0][1] = 1.0;    m[0][2] = 0.0;
  1098.     m[n][0] = 0.0;    m[n][1] = 1.0;    m[n][2] = 0.0;
  1099.     r[0] = ty1;                           /* Set start tangent. */
  1100.     r[n] = ty2;                         /* Set end tangent. */
  1101.     t = t_min;
  1102.     dt = (t_max-t_min)/n;
  1103.     for (i=1; i<n; i++) {
  1104.        t = t + dt;
  1105.        m[i][0] = dt;
  1106.        m[i][2] = dt;
  1107.        m[i][1] = 2 * (m[i][0] + m[i][2]);
  1108.        r[i] = m[i][0] * ((p -> next -> Y) - (p -> Y)) / m[i][2]
  1109.             + m[i][2] * ((p -> next -> next -> Y) -
  1110.                          (p -> next -> Y)) / m[i][0];
  1111.        r[i] *= 3.0;
  1112.        p = p -> next;
  1113.     }
  1114.  
  1115.     if (!solve_tri_diag(m, r, tangents_y, n+1)) { /* Find the Y(t) tangents. */
  1116.         free((char *) m);
  1117.         free((char *) r);
  1118.     int_error("Cannt interpolate Y using complete splines", NO_CARET);
  1119.     }
  1120.     free((char *) m);
  1121.     free((char *) r);
  1122. }
  1123.  
  1124. /*
  1125.  * Solve tri diagonal linear system equation. The tri diagonal matrix is
  1126.  * defined via matrix M, right side is r, and solution X i.e. M * X = R.
  1127.  * Size of system given in n. Return TRUE if solution exist.
  1128.  */
  1129. int solve_tri_diag(m, r, x, n)
  1130. tri_diag m[];
  1131. double r[], x[];
  1132. int n;
  1133. {
  1134.     int i;
  1135.     double t;
  1136.  
  1137.     for (i=1; i<n; i++) {   /* Eliminate element m[i][i-1] (lower diagonal). */
  1138.     if (m[i-1][1] == 0) return FALSE;
  1139.     t = m[i][0] / m[i-1][1];        /* Find ratio between the two lines. */
  1140. /*    m[i][0] = m[i][0] - m[i-1][1] * t; */
  1141. /* m[i][0] is not used any more (and set to 0 in the above line) */
  1142.     m[i][1] = m[i][1] - m[i-1][2] * t;
  1143.     r[i] = r[i] - r[i-1] * t;
  1144.     }
  1145.     /* Now do back subtitution - update the solution vector X: */
  1146.     if (m[n-1][1] == 0) return FALSE;
  1147.     x[n-1] = r[n-1] / m[n-1][1];               /* Find last element. */
  1148.     for (i=n-2; i>=0; i--) {
  1149.     if (m[i][1] == 0) return FALSE;
  1150.     x[i] = (r[i] - x[i+1] * m[i][2]) / m[i][1];
  1151.     }
  1152.     return TRUE;
  1153. }
  1154.  
  1155. /*
  1156.  * Generate a Bspline curve defined by all the points given in linked list p:
  1157.  * Algorithm: using deBoor algorithm
  1158.  * Note: if Curvekind is OPEN_CONTOUR than Open end knot vector is assumed,
  1159.  *       else (CLOSED_CONTOUR) Float end knot vector is assumed.
  1160.  * It is assumed that num_of_points is at list 2, and order of Bspline is less
  1161.  * than num_of_points!
  1162.  */
  1163. static void gen_bspline_approx(p_cntr, num_of_points, order, contour_kind)
  1164. struct cntr_struct *p_cntr;
  1165. int num_of_points, order, contour_kind;
  1166. {
  1167.     int i, knot_index = 0, pts_count = 1;
  1168.     double dt, t, next_t, t_min, t_max, x, y;
  1169.     struct cntr_struct *pc_temp = p_cntr, *pc_tail=NULL;
  1170.  
  1171.     /* If the contour is Closed one we must update few things:               */
  1172.     /* 1. Make the list temporary circular, so we can close the contour.     */
  1173.     /* 2. Update num_of_points - increase it by "order-1" so contour will be */
  1174.     /*    closed. This will evaluate order more sections to close it!        */
  1175.     if (contour_kind == CLOSED_CONTOUR) {
  1176.     pc_tail = p_cntr;
  1177.     while (pc_tail -> next) pc_tail = pc_tail -> next;/* Find last point.*/
  1178.     pc_tail -> next = p_cntr;   /* Close contour list - make it circular.*/
  1179.     num_of_points += order;
  1180.     }
  1181.  
  1182.     /* Find first (t_min) and last (t_max) t value to eval: */
  1183.     t = t_min = fetch_knot(contour_kind, num_of_points, order, order);
  1184.     t_max = fetch_knot(contour_kind, num_of_points, order, num_of_points);
  1185.     next_t = t_min + 1.0;
  1186.     knot_index = order;
  1187.     dt = 1.0/num_approx_pts;            /* Number of points per one section. */
  1188.  
  1189.  
  1190.     while (t<t_max) {
  1191.     if (t > next_t) {
  1192.         pc_temp = pc_temp -> next;     /* Next order ctrl. pt. to blend. */
  1193.             knot_index++;
  1194.         next_t += 1.0;
  1195.     }
  1196.         eval_bspline(t, pc_temp, num_of_points, order, knot_index,
  1197.                         contour_kind, &x, &y);   /* Next pt. */
  1198.     add_cntr_point(x, y);
  1199.     pts_count++;
  1200.     /* As we might have some real number round off problems we must      */
  1201.     /* test if we dont produce too many points here...                   */
  1202.     if (pts_count + 1 == num_approx_pts * (num_of_points - order) + 1)
  1203.             break;
  1204.         t += dt;
  1205.     }
  1206.  
  1207.     eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index,
  1208.         contour_kind, &x, &y);
  1209.     /* If from round off errors we need more than one last point: */
  1210.     for (i=pts_count; i<num_approx_pts * (num_of_points - order) + 1; i++)
  1211.     add_cntr_point(x, y);                /* Complete the contour. */
  1212.  
  1213.     if (contour_kind == CLOSED_CONTOUR)     /* Update list - un-circular it. */
  1214.     pc_tail -> next = NULL;
  1215. }
  1216.  
  1217. /*
  1218.  * The recursive routine to evaluate the B-spline value at point t using
  1219.  * knot vector PKList, and the control points Pdtemp. Returns x, y after the
  1220.  * division by the weight w. Note that Pdtemp points on the first control
  1221.  * point to blend with. The B-spline is of order order.
  1222.  */
  1223. static void eval_bspline(t, p_cntr, num_of_points, order, j, contour_kind, x, y)
  1224. double t;
  1225. struct cntr_struct *p_cntr;
  1226. int num_of_points, order, j, contour_kind;
  1227. double *x, *y;
  1228. {
  1229.     int i, p;
  1230.     double ti, tikp, *dx, *dy;      /* Copy p_cntr into it to make it faster. */
  1231.  
  1232.     dx = (double *) alloc((unsigned long) (sizeof(double) * (order + j)),
  1233.                         "contour b_spline");
  1234.     dy = (double *) alloc((unsigned long) (sizeof(double) * (order + j)),
  1235.                         "contour b_spline");
  1236.     /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */
  1237.     for (i=j-order; i<=j; i++) {
  1238.         dx[i] = p_cntr -> X;
  1239.         dy[i] = p_cntr -> Y;
  1240.         p_cntr = p_cntr -> next;
  1241.     }
  1242.  
  1243.     for (p=1; p<=order; p++) {        /* Iteration (b-spline level) counter. */
  1244.     for (i=j; i>=j-order+p; i--) {           /* Control points indexing. */
  1245.             ti = fetch_knot(contour_kind, num_of_points, order, i);
  1246.             tikp = fetch_knot(contour_kind, num_of_points, order, i+order+1-p);
  1247.         if (ti == tikp) {   /* Should not be a problems but how knows... */
  1248.         }
  1249.         else {
  1250.         dx[i] = dx[i] * (t - ti)/(tikp-ti) +         /* Calculate x. */
  1251.             dx[i-1] * (tikp-t)/(tikp-ti);
  1252.         dy[i] = dy[i] * (t - ti)/(tikp-ti) +         /* Calculate y. */
  1253.             dy[i-1] * (tikp-t)/(tikp-ti);
  1254.         }
  1255.     }
  1256.     }
  1257.     *x = dx[j]; *y = dy[j];
  1258.     free((char *) dx);
  1259.     free((char *) dy);
  1260. }
  1261.  
  1262. /*
  1263.  * Routine to get the i knot from uniform knot vector. The knot vector
  1264.  * might be float (Knot(i) = i) or open (where the first and last "order"
  1265.  * knots are equal). contour_kind determines knot kind - OPEN_CONTOUR means
  1266.  * open knot vector, and CLOSED_CONTOUR selects float knot vector.
  1267.  * Note the knot vector is not exist and this routine simulates it existance
  1268.  * Also note the indexes for the knot vector starts from 0.
  1269.  */
  1270. static double fetch_knot(contour_kind, num_of_points, order, i)
  1271. int contour_kind, num_of_points, order, i;
  1272. {
  1273.     switch (contour_kind) {
  1274.     case OPEN_CONTOUR:
  1275.         if (i <= order) return 0.0;
  1276.         else if (i <= num_of_points) return (double) (i - order);
  1277.          else return (double) (num_of_points - order);
  1278.     case CLOSED_CONTOUR:
  1279.         return (double) i;
  1280.     default: /* Should never happen */
  1281.         return 1.0;
  1282.     }
  1283. #ifdef sequent
  1284.     return 1.0;
  1285. #endif
  1286. }
  1287.